# QJPS Replication
# Figures 1 and 2
# Updated by David K. Park

# Set Working Directory HERE
setwd("C:/Users/F001CKB/Documents/Political Effects of Land-Use Regulation/gelman data/QJPS_6026_supp/")

# load foreign package to read in STATA file
library(foreign)

fips.cbs <- read.dta("fips.icpsr.cbs.naes.dta")

# Figure 1 Regressions
# state-level
st <- read.csv("st.csv")

st.all <- merge(st, fips.cbs, by.x="stateabb", by.y="st")
st.s <- subset(st.all, st.all$regioncbs==3)
st.ns <- subset(st.all, st.all$regioncbs!=3)

p.st.yr <- seq(min(st$year), max(st$year),4)
n.st.yr <- length(p.st.yr)

col.st.names <- c("st.int", "stcpi")
  n.st.coef <- length(col.st.names)

B.st.lm <- array(NA, c(n.st.yr, n.st.coef))
  colnames(B.st.lm) <- col.st.names
SE.st.lm <- array(NA, c(n.st.yr, n.st.coef))
  colnames(SE.st.lm) <- col.st.names

B.st.s.lm <- array(NA, c(n.st.yr, n.st.coef))
  colnames(B.st.s.lm) <- col.st.names
SE.st.s.lm <- array(NA, c(n.st.yr, n.st.coef))
  colnames(SE.st.s.lm) <- col.st.names

B.st.ns.lm <- array(NA, c(n.st.yr, n.st.coef))
  colnames(B.st.ns.lm) <- col.st.names
SE.st.ns.lm <- array(NA, c(n.st.yr, n.st.coef))
  colnames(SE.st.ns.lm) <- col.st.names

for (i in 1:n.st.yr){
    st.temp <- subset(st, st$year==p.st.yr[i])
    st.reg <- summary(lm(st.temp$st.repshare ~ st.temp$stcpi))
    B.st.lm[i,] <- st.reg$coef[,1]
    SE.st.lm[i,] <- st.reg$coef[,2]
}

for (i in 1:n.st.yr){
    st.s.temp <- subset(st.s, st.s$year==p.st.yr[i])
    st.s.reg <- summary(lm(st.s.temp$st.repshare ~ st.s.temp$stcpi))
    B.st.s.lm[i,] <- st.s.reg$coef[,1]
    SE.st.s.lm[i,] <- st.s.reg$coef[,2]
}

for (i in 1:n.st.yr){
    st.ns.temp <- subset(st.ns, st.ns$year==p.st.yr[i])
    st.ns.reg <- summary(lm(st.ns.temp$st.repshare ~ st.ns.temp$stcpi))
    B.st.ns.lm[i,] <- st.ns.reg$coef[,1]
    SE.st.ns.lm[i,] <- st.ns.reg$coef[,2]
}

B.st.lm <- data.frame(B.st.lm)
SE.st.lm <- data.frame(SE.st.lm)

B.st.s.lm <- data.frame(B.st.s.lm)
SE.st.s.lm <- data.frame(SE.st.s.lm)

B.st.ns.lm <- data.frame(B.st.ns.lm)
SE.st.ns.lm <- data.frame(SE.st.ns.lm)

#  Figure 1 Plot
windows(width=8, height=2.5)
par(mfrow=c(1,3))
n.st <- NROW(B.st.lm)

plot(0, 0, type='n', ylim=c(-.5, .5), xlim=c(min(st$year), max(st$year)),
  cex.lab=1.15, xaxt="n", yaxt="n",  main="All States", cex.main=1, xlab="Year", ylab="CPI Coefficient")
  axis(side=1, at=c(1960, 1980, 2000), cex.axis=1.1)
  axis(side=2, at=c(-.4, -.2, 0, .2, .4), cex.axis=1.1)
  abline(h=0, col="gray")
    for (i in 1:n.st) {
    points(p.st.yr[i], B.st.lm[i,2], type='p', pch=19)
    segments(p.st.yr[i], (B.st.lm[i,2]-SE.st.lm[i,2]), p.st.yr[i], (B.st.lm[i,2]+SE.st.lm[i,2]))
  }

plot(0, 0, type='n', ylim=c(-.75, .5), xlim=c(min(st$year), max(st$year)),
  cex.lab=1.15, xaxt="n", yaxt="n",  main="Southern States", cex.main=1, xlab="Year", ylab="CPI Coefficient")
  axis(side=1, at=c(1960, 1980, 2000), cex.axis=1.1)
  axis(side=2, at=c(-.6, -.4, -.2, 0, .2, .4), cex.axis=1.1)
  abline(h=0, col="gray")
    for (i in 1:n.st) {
    points(p.st.yr[i], B.st.s.lm[i,2], type='p', pch=19)
    segments(p.st.yr[i], (B.st.s.lm[i,2]-SE.st.s.lm[i,2]), p.st.yr[i], (B.st.s.lm[i,2]+SE.st.s.lm[i,2]))
  }

plot(0, 0, type='n', ylim=c(-.5, .5), xlim=c(min(st$year), max(st$year)),
  cex.lab=1.15, xaxt="n", yaxt="n",  main="Non-Southern States", cex.main=1, xlab="Year", ylab="CPI Coefficient")
  axis(side=1, at=c(1960, 1980, 2000), cex.axis=1.1)
  axis(side=2, at=c(-.4, -.2, 0, .2, .4), cex.axis=1.1)
  abline(h=0, col="gray")
    for (i in 1:n.st) {
    points(p.st.yr[i], B.st.ns.lm[i,2], type='p', pch=19)
    segments(p.st.yr[i], (B.st.ns.lm[i,2]-SE.st.ns.lm[i,2]), p.st.yr[i], (B.st.ns.lm[i,2]+SE.st.ns.lm[i,2]))
  }

# Figure 1
savePlot(filename="st.yr", type=c("pdf"), device=dev.cur())
savePlot(filename="st.yr", type=c("ps"), device=dev.cur())

